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By making use of the data on the total signal and on the muon component of the air showers detected 
by the Yakutsk array, we analyze, in the frameworks of the recently suggested event-by-event approach, how 
large the fraction of primary gamma-rays at ultra-high energies can be. We derive upper limits on the photon 
fraction in the integral flux of primary cosmic rays. At the 95% confidence level (CL), these limits are 22% 
for primary energies E > 4 ■ 10 19 eV and 12% for E > 2 ■ 10 19 eV. Despite the presence of muonless events, 
the data are consistent with complete absence of photons at least at 95% CL. The sensitivity of the results 
to systematic uncertainties, in particular to those of the energy determination for non-photon primaries, is 
discussed. 
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With the increase of statistics of ultra-high-energy 
(UHE) cosmic-ray (CR) events the study of the chemi- 
cal composition at the very end of the spectrum (beyond 
10 19 eV) is becoming quite realistic. This issue is of 
primary interest today, in view of a systematic discrep- 
ancy between energy spectra measured by different de- 
tectors [UmElll] an d °f indications towards a fraction of 
>- neutral particles among the UHECR primaries [5] . The 
chemical composition is also a starting point in stud- 
5_j des of: (i) extragalactic magnetic fields and radiation 
^backgrounds; (ii) accelaration mechanisms operating in 
astrophysical sources; (in) possible top-down scenarios 
emerging in various extensions of the Standard Model 
of particle physics. In particular, the photon fraction 
in the CR flux is of crucial importance; the aim of this 
work is to derive stringent limits on this fraction in the 
integral CR flux above the energy 2 • 10 19 eV. 

We make use of a recently suggested approach [51 [7] 
and perform case-by-case analysis of 50 events detected 
by the Yakutsk extensive-air-shower array (Yakutsk ar- 
ray in what follows) lj with reconstructed energies 
above 2 • 10 19 eV chosen according to quality cuts de- 
scribed in Sec. [2j To place the limit on the photon 
fraction, we compare the reported information on sig- 
nals measured by scintillation and muon detectors with 
that expected from air-shower simulations. We focus on 
the surface detector signal density at 600 meters 5(600) 
and the muon density at 1000 meters, p M (1000), which 
are used in experiments as primary energy and primary 



composition estimators, respectively. Among the fifty 
showers in the sample, two are muonless (that is, muon 
detectors were operating in the shower impact area but 
did nor detect any signal) . These events are compatible 
with being initiated by primary gamma rays of energies 
2 • 10 19 eV< E Q < 4 ■ 10 19 eV, even though the recon- 
structed energy exceeds 4 • 10 19 eV for one of them. One 
muon-poor shower is consistent with a photon primary 
of energy above 4 • 10 19 eV with probability about 10%. 
For the rest of the showers, the hypothesis of a photon 
primary is rejected at the 95% CL for each event. We 
derive upper limits on the fraction e 7 of photons in the 
integral flux of primary cosmic rays with actual energies 
E Q > 2 • 10 19 eV and E > 4 • 10 19 eV (the difference 
between actual (En) and reconstructed (E es t) energies 
is discussed in Sec. [3]). 

The rest of the paper is organized as follows. In 
Sec. [2] we discuss the experimental data set used in our 
study. In Sec. [3] we briefly review the approach we use 
and present our main results. We discuss how robust 
these results are with respect to changes in assumptions 
and in the analysis procedure and discuss the uncertain- 
ties associated with possible systematics in energy de- 
termination of observed UHECR events in Sec.|H Sec. [5] 
contains our conclusions. 



2 EXPERIMENTAL DATA 

Yakutsk array is observing UHECR events since 
1973, with detectors in various configurations. Since 
1979, muon detectors with areas up to 36 m 2 (cur- 



1 



2 



A.V. Glushkov et al. 



rently, five detectors of 20 m 2 each with threshold energy 
1 GeV for vertical muons) supplement ground-based 
scintillator stations. At present, it is the only installa- 
tion equipped with muon detectors capable of studying 
ultra-high-energy cosmic rays. 

The energy of a primary particle is estimated from 
5(600) and zenith angle with the help of the procedure 
described in Ref. [8] , calibrated experimentally by mak- 
ing use of the atmospheric Cherenkov light. This re- 
constructed energy E cst differs from the true primary 
energy Eq both due to natural fluctuations and due to 
possible systematic effects. These latter effects depend 
on the primary particle type; in particular, the differ- 
ence between photons and hadrons is significant. More- 
over, for photons, the effects of geomagnetic field [9] re- 
sult in directional dependence of the energy reconstruc- 
tion. Thus, the event energy reported by the experiment 
should be treated with care when we allow the primary 
to be a photon. Because of possible energy underesti- 
mation for high-energy photon-induced showers, we use 
events with E cst > 2 x f 19 eV even when deriving the 
limit for Eq > 4 • 10 19 eV; they contribute to the final 
limit with different weights [7]. 

For our study, we selected a subset of events with 
E cs t > 2 x fO 19 eV satisfying the following cuts aimed 
at the most precise determination of both 5(600) and 
M1000): 

(i) shower core inside the array; 

(it) zenith angle 9 < 60°; 

(Hi) three or more muon detectors between 400 m 
and 2000 m from the shower axis, operational at the 
moment of the shower arrival. 

Our sample consists of 50 air showers; the cuts se- 
lect approximately one third of the events used for the 
determination of the spectrum. 



3 SIMULATIONS AND RESULTS 

The approach we use was described and discussed in 
detail in Refs. [6] [7] and has already been applied to a 
similar study of the photon fraction at energies above 
fO 20 eV [6]. Here, we summarise the main steps of this 
approach. 

For each of the events in our sample, we generated 
a library of simulated showers induced by primary pho- 
tons. Thrown energies Eq of the simulated showers were 
randomly selected within a relevant energy interval in 
order to take into account possible deviations of E est 
from Eq, see below. The arrival directions of the simu- 
lated showers were the same as those of the correspond- 
ing real events. The simulations were performed with 



CORS1KA v6.5011 [TU], choosing QGSJET 11-03 dJ 
as high-energy and FLUKA 2005.6 12\ as low-energy 
hadronic interaction model. Electromagnetic shower- 
ing was implemented with EGS4 [13] incorporated into 
CORSIKA. Possible interactions of the primary pho- 
tons with the geomagnetic field were simulated with the 
PRESHOWER option of CORSIKA 14J. As suggested 
in Ref. [15j . all simulations were performed with thin- 
ning level 10~ 5 , maximal weight 10 6 for electrons and 
photons, and fO 4 for hadrons. 

For each simulated shower, we determined 5(600) 
and /^(fOOO) by making use of the detector response 
functions from Ref. [16]. For a given arrival direction, 
there is one-to-one correspondence between 5(600) and 
the estimated energy as determined by the standard 
analysis procedure for the Yakutsk experiment [8] . This 
enables us to select simulated showers compatible with 
the observed ones by the signal density, which follows 
the Gaussian distribution in log(energy); the standard 
deviation of E est has been determined event-by-event 
and is typically 17% [17] • Namely, to each simulated 
shower, we assigned a weight w\ proportional to this 
Gaussian probability distribution in log E cs t centered at 
the observed energy. Additional weight u>2 was assigned 
to each simulated shower to reproduce the thrown en- 
ergy spectrum oc E^ 2 (see Sec. 14. 31 for the discussion of 
the variation of the spectral index). For each of the ob- 
served events from our dataset, we calculated the distri- 
bution of muon densities p M (1000) representing photon- 
induced showers compatible with the observed ones by 
5(600) and arrival directions. To this end, we calculated 
p M (1000) for each simulated shower by making use of the 
same muon lateral distribution function as used in the 
analysis of real data [18] . To take into account possible 
experimental errors in the determination of the muon 
density, we replaced each simulated p M (1000) by a distri- 
bution representing possible statistical errors (Gaussian 
with 25% standard deviation [5]). The distribution of 
the simulated muon densities is the sum of these Gaus- 
sians weighted by w\W2- 

For each event we calculate, by making use of the ob- 
tained distributions, two numbers: the probability that 
it could be initiated by a photon with true energy in 
the range of interest (that is, above Eq = 4 ■ 10 19 eV 
or above Eq — 2 • 10 19 eV) and the probability that it 
could be initiated by any other primary (whose energy 
is assumed to be determined correctly by the experi- 
ment; see Sec. 14.11 for relaxing this assumption) with 
energy above this Eq. For most of the events, the mea- 
sured muon densities are too high as compared to those 
obtained from simulations of photon induced showers. 

Given these probabilities for each event, we con- 
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Figure 1. Limits (95% CL) on the fraction e 7 of photons 
in the integral CR flux versus energy. The results of the 
present work (Y) are shown together with the limits pre- 
viously given in Refs. [19] (HP), [20] (A), [6] (AY) and 
[21] (PA). Also shown are predictions for the superheavy 
dark matter model (thick line), the topological-defect 
models (necklaces, between dotted lines) |22| and the 
Z-burst model (shaded area) 23 . Theoretical curves 
are normalized to the AGASA spectrum [2]. Energy is 
measured in eV. 

struct the likelihood function (see Ref. [7J for details) 
to estimate, at a given confidence level, the fraction e 7 
of primary photons among UHECR with energies in a 
given range. In this way we obtain at 95% CL 

e 7 < 22% for E > 4 • 10 19 eV, (1) 

e 7 < 12% for E > 2 • 10 19 eV. (2) 

These limits include corrections for the "lost photons" 
(those with true energies Eq > 4-10 19 eV for the limit |T|) 
and Eq > 2 • 10 19 eV for the limit (|2|) but reconstructed 
energies E YCC < 2 • 10 19 eV, see Ref. [TJ for more details). 

In Fig. [TJ we present our limits (denoted by Y) to- 
gether with previously published limits 1 ) on e 7 . Also, 
typical theoretical predictions are shown for the super- 
heavy dark matter, topological-defect and Z-burst mod- 
els. Our limits on e 7 are currently the strongest ones for 
the energy range under discussion. They disfavor the su- 
perheavy dark matter explanation of the highest energy 
events. 



4 ROBUSTNESS OF THE RESULTS 

The systematic uncertainties of our results are re- 
lated to the air-shower simulations and to the data in- 
terpretation. They were discussed in detail in Ref. [6] 
for a different data set, with the conclusion that the 
approach we use to constrain e 7 results in quite robust 
limits. 

4.1 Systematic uncertainty in the 5(600) and 
energy determination 

The systematic uncertainty in the absolute energy 
determination by the Yakutsk array is about 30% [JJ. 
It originates from two quite different sources: (a) the 
measurement of 5(600) and (b) the relation between 
5(600) and primary energy. The probabilities that a 
particular event may allow for a gamma-ray interpreta- 
tion are not at all sensitive to the 5(600)-to-energy con- 
version because we select simulated events by 5(600) 
and not by energy. These probabilities may only be 
affected by relative systematics between the determina- 
tions of /o^(1000) and of 5(600). On the other hand, we 
assumed that the experimental energy determination is 
correct for non-photon primaries; the values of proba- 
bilities that a particular event could be initiated by a 
non-photon primary with energy above threshold and 
hence the effective number of events contributing to the 
limit on e 7 would change if the energies are systemat- 
ically shifted. The effect of such a shift would be to 
change the energy range for which the limit is applica- 
ble and to change, by a few per cent, the limit itself. 
This is illustrated in Fig. [21 which has been obtained 
in a way similar to that described in Sec. El but with 
six minimal values of Eq for each of the three curves 
corresponding to —30%, 0% and +30% shifts in ener- 
gies of non-photon primaries. We see that the limit at 
Eq > 2 ■ 10 19 eV is uncertain by less than a few per cent 
while at higher energies, systematic shifts downwards 
reduce statistics considerably, which results in relaxing 
the limit. Similar uncertainties are expected for limits 
from other experiments shown in Fig. [TJ Note that the 
theoretical expectations presented there are also sensi- 
tive to the energy scale. 



!'A 65% upper limit for energies above 1.2 ■ 10 20 eV has been 
claimed from the study of AGASA data [24] : however, there are 
problems in accounting for the difference between actual and re- 
constructed photon energies in that work (see Ref. [6] for a de- 
tailed discussion). 



4.2 Interaction models and simulation codes 

Our simulations were performed entirely in the 
CORSIKA framework, and any change in the interac- 
tion models or simulation codes, which affects either 
5(600) or /^(lOOO), may affect our limit. As discussed 
in Ref. 6J, our method is quite robust with respect to 
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Figure 2. Sensitivity of the limit on e 7 to the systematic 
uncertainty in the energy determination for non-photon 
primaries. The solid (red) curve represents the limits 
assuming the energy scale quoted by Yakutsk experi- 
ment and normalized to the Cherenkov light (the same 
as points Y in FigfJ}. The dashed (blue) and dotted 
(black) curves correspond to the shifts of —30% and 
+30%, respectively, in all energies of non- photon pri- 
maries. 

the changes in the interaction models and to reason- 
able variations in the extrapolation of the photonuclear 
cross section to high energies (the values presented here 
were obtained for the standard parameterization of the 
photonuclear cross section given by the Particle Data 
Group [35] and implemented as default in CORSIKA). 

4.3 Primary energy spectrum 

For our limit, we used the primary photon spectrum 
Eq" for a = 2. Change in the value of a affects the 
final limit on e 7 through the fraction of "lost" photons, 
but we have found that variations of a in the interval 
1 < a < 3 result in variations of e 7 only within 1%. 

4.4 Width of the distribution 

The rare probabilities of high values of /^(lOOO) in 
the tail of the distribution for primary photons depend 
on the width of this distribution. The following sources 
contribute to this width: 

• variations of the primary energy compatible with 
the observed 5(600) (larger energy corresponds to 
larger muon number and hence /3 M (I000)); 

• physical shower-to-shower fluctuations in muon 
density for a given energy (dominated by fluc- 
tuations in the first few interactions, including 
preshowering in the geomagnetic field); 



• artificial fluctuations in 5(600) and ^(1000) due 
to thinning; 

• experimental errors in /^(fOOO) determination. 

While the first two sources are physical and are fully 
controlled by the simulation code, the variations of the 
last two may affect the results. 

It has been noted in Ref. [53] that the fluctuations in 
/9 M (I000) due to thinning may affect strongly the preci- 
sion of the composition studies. For the thinning param- 
eters we use, the relative size of these fluctuations [UJ 
is - 10% for /v(I000) and ~ 5% for 5(600). Thus with 
more precise simulations, the distributions of muon den- 
sities should become more narrow, which would reduce 
the probability of the gamma-ray interpretation of the 
studied events even further. 

The distributions of ,0^(1000) we use accounted for 
the error in the experimental determination of this 
quantity. In principle, this error depends on the event 
quality and on the muon number itself, which is sys- 
tematically lower for simulated gamma-induced show- 
ers than for the observed events. Still, we tested the 
stability of our limit by taking artificially high values 
of experimental errors in muon density: 50% instead of 
25%. The limit on e 7 changes by less than one per cent 
in that case. 



5 CONCLUSIONS 

To summarize, we have studied the possibility that 
ultra-high energy events observed by the Yakutsk array 
were initiated by primary photons. The use of large-area 
muon detectors, a unique feature of the Yakutsk exper- 
iment, together with the new analysis method [6] [7], 
enabled us to put stringent constraints on the gamma- 
ray primaries even with a relatively small set of high- 
quality data. An important ingredient in our study was 
the careful tracking of differences between the actual 
and reconstructed energies. We obtained upper bounds 
([!]), ([2]) on the fraction e 7 of primary photons, assum- 
ing an isotropic photon flux and Eq 2 spectrum. These 
limits are the strongest ones up to date; they constrain 
considerably the superheavy dark matter models. 
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